function [best_schedule, optimization_history, total_time] = grid_saea_rfs(grid_data, algorithm_params)
% 基于SAEA-RFS算法的智能电网调度优化
% 将电网调度问题建模为高维优化问题，使用SAEA-RFS算法求解
%
% INPUTS:
%   grid_data - 电网模型数据
%   algorithm_params - 算法参数
%
% OUTPUTS:
%   best_schedule - 最优调度方案
%   optimization_history - 优化历史记录
%   total_time - 总运行时间

% 屏蔽RBF矩阵奇异值警告
warning('off', 'MATLAB:nearlySingularMatrix');
warning('off', 'MATLAB:singularMatrix');

%% 算法参数设置
if nargin < 2
    algorithm_params = struct();
end

% 默认参数
if ~isfield(algorithm_params, 'population_size'), algorithm_params.population_size = 20; end
if ~isfield(algorithm_params, 'initial_samples'), algorithm_params.initial_samples = 100; end
if ~isfield(algorithm_params, 'sub_problems'), algorithm_params.sub_problems = 15; end
if ~isfield(algorithm_params, 'sub_problem_size'), algorithm_params.sub_problem_size = 80; end
if ~isfield(algorithm_params, 'max_iterations'), algorithm_params.max_iterations = 10; end
if ~isfield(algorithm_params, 'max_evaluations'), algorithm_params.max_evaluations = 5000; end

%% 问题维度计算
num_gen = grid_data.num_generators;
num_storage = grid_data.num_storage;
num_loads = grid_data.num_loads;
time_periods = grid_data.time_periods;

% 决策变量维度
D = num_gen * time_periods + num_storage * time_periods + num_loads * time_periods;
fprintf('电网调度问题维度: %d\n', D);
fprintf('- 发电机出力变量: %d\n', num_gen * time_periods);
fprintf('- 储能充放电变量: %d\n', num_storage * time_periods);
fprintf('- 负荷分配变量: %d\n', num_loads * time_periods);

%% 搜索空间边界设置
% 发电机出力边界
gen_bounds = [];
for g = 1:num_gen
    gen = grid_data.generators(g);
    gen_bounds = [gen_bounds, repmat([gen.min_power, gen.max_power], 1, time_periods)];
end

% 储能设备功率边界
storage_bounds = [];
for s = 1:num_storage
    storage = grid_data.storage(s);
    storage_bounds = [storage_bounds, repmat([-storage.max_power, storage.max_power], 1, time_periods)];
end

% 负荷分配边界（0.8-1.2倍基础负荷）
load_bounds = [];
for l = 1:num_loads
    load_bounds = [load_bounds, repmat([0.8, 1.2], 1, time_periods)];
end

% 合并所有边界
lu = [gen_bounds, storage_bounds, load_bounds; ...
      gen_bounds, storage_bounds, load_bounds];
lu = reshape(lu, 2, D);

%% 初始化训练样本档案
fprintf('初始化训练样本...\n');
archive = zeros(algorithm_params.initial_samples, D + 1);

% 使用拉丁超立方采样生成初始解
for i = 1:algorithm_params.initial_samples
    % 生成随机解
    x = lu(1, :) + (lu(2, :) - lu(1, :)) .* rand(1, D);
    
    % 评估目标函数
    [cost, ~] = grid_objective(x, grid_data);
    
    % 存储到档案
    archive(i, 1:D) = x;
    archive(i, D+1) = cost;
end

% 找到初始最优解
[best_cost, best_idx] = min(archive(:, D+1));
best_solution = archive(best_idx, 1:D);

fprintf('初始最优成本: %.2f\n', best_cost);

%% 初始化种群
fprintf('初始化优化种群...\n');
population = zeros(algorithm_params.population_size, D);
for i = 1:algorithm_params.population_size
    population(i, :) = lu(1, :) + (lu(2, :) - lu(1, :)) .* rand(1, D);
end

%% 主优化循环
fprintf('开始SAEA-RFS优化...\n');
evaluation_count = algorithm_params.initial_samples;
optimization_history = [evaluation_count, best_cost];
cycle = 0;

tic;
while evaluation_count < algorithm_params.max_evaluations
    cycle = cycle + 1;
    fprintf('优化周期 %d: 当前最优成本 = %.2f, 评估次数 = %d\n', ...
        cycle, best_cost, evaluation_count);
    
    %% 1. 使用RFS构建子问题并训练代理模型
    [lambda, gamma, dv, sp] = rfs(algorithm_params.sub_problem_size, ...
        archive(:, 1:D), archive(:, D+1), 'cubic', algorithm_params.sub_problems);
    
    %% 2. 子问题优化
    sub_eval_min = [];
    
    for i = 1:algorithm_params.sub_problems
        % 提取子问题种群
        sub_population = population(:, dv{i});
        
        % 利用RBF模型预测子问题种群适应度
        sub_fitness = rbf_eval(sub_population, archive(sp{i}, dv{i}), ...
            lambda{i}, gamma{i}, 'cubic');
        
        % 子问题搜索空间边界
        sub_bounds = lu(:, dv{i});
        
        % 使用DE优化该子问题
        for iter = 1:algorithm_params.max_iterations
            % DE生成新候选解
            sub_offspring = de(sub_population, best_solution(dv{i}), sub_bounds);
            
            % 利用RBF模型预测子代适应度
            sub_offspring_fitness = rbf_eval(sub_offspring, archive(sp{i}, dv{i}), ...
                lambda{i}, gamma{i}, 'cubic');
            
            % 合并父代与子代并排序
            all_fitness = [sub_fitness; sub_offspring_fitness];
            all_population = [sub_population; sub_offspring];
            
            [~, sorted_idx] = sort(all_fitness);
            
            % 选择最优的个体作为新的子问题种群
            sub_population = all_population(sorted_idx(1:algorithm_params.population_size), :);
            sub_fitness = all_fitness(sorted_idx(1:algorithm_params.population_size));
        end
        
        % 更新主种群
        population(:, dv{i}) = sub_population;
        
        % 记录子问题最优解
        [min_fitness, min_idx] = min(sub_fitness);
        sub_eval_min = [sub_eval_min; min_fitness, min_idx, i];
    end
    
    %% 3. 真实函数评估
    % 选择最有希望的候选解进行真实评估
    num_real_eval = min(10, algorithm_params.max_evaluations - evaluation_count);
    
    % 从所有子问题的最优解中选择
    [~, top_idx] = sort(sub_eval_min(:, 1));
    top_sub_problems = sub_eval_min(top_idx(1:num_real_eval), 3);
    
    for i = 1:num_real_eval
        sub_idx = top_sub_problems(i);
        candidate_idx = sub_eval_min(sub_eval_min(:, 3) == sub_idx, 2);
        
        % 构建完整解
        candidate_solution = best_solution; % 从当前最优解开始
        candidate_solution(dv{sub_idx}) = population(candidate_idx, dv{sub_idx});
        
        % 真实函数评估
        [real_cost, ~] = grid_objective(candidate_solution, grid_data);
        evaluation_count = evaluation_count + 1;
        
        % 更新档案
        archive = [archive; candidate_solution, real_cost];
        
        % 更新全局最优解
        if real_cost < best_cost
            best_cost = real_cost;
            best_solution = candidate_solution;
            fprintf('发现新的最优解: %.2f (评估次数: %d)\n', best_cost, evaluation_count);
        end
        
        % 记录优化历史
        optimization_history = [optimization_history; evaluation_count, best_cost];
        
        % 检查是否达到最大评估次数
        if evaluation_count >= algorithm_params.max_evaluations
            break;
        end
    end
    
    % 限制档案大小
    if size(archive, 1) > 1000
        [~, sorted_idx] = sort(archive(:, D+1));
        archive = archive(sorted_idx(1:1000), :);
    end
end

total_time = toc;

%% 结果整理
best_schedule = struct();
best_schedule.generator_power = reshape(best_solution(1:num_gen*time_periods), num_gen, time_periods);
best_schedule.storage_power = reshape(best_solution(num_gen*time_periods+1:num_gen*time_periods+num_storage*time_periods), num_storage, time_periods);
best_schedule.load_allocation = reshape(best_solution(num_gen*time_periods+num_storage*time_periods+1:end), num_loads, time_periods);
best_schedule.total_cost = best_cost;

fprintf('\n优化完成！\n');
fprintf('总运行时间: %.2f 秒\n', total_time);
fprintf('总评估次数: %d\n', evaluation_count);
fprintf('最终最优成本: %.2f\n', best_cost);

% 恢复警告设置
warning('on', 'MATLAB:nearlySingularMatrix');
warning('on', 'MATLAB:singularMatrix');

end

